perm filename ARITH.PAL[2,VDS] blob sn#163104 filedate 1975-12-29 generic text, type C, neo UTF8
COMMENT ⊗   VALID 00012 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00003 00002	.TITLE ARITH
C00005 00003	ASIN - FLOATING POINT, SINGLE PRECISION ARC-SINE FUNCTION
C00007 00004	ACOS - FLOATING POINT, SINGLE PRECISION ARC-COSINE FUNCTION
C00009 00005	ATAN - FLOATING POINT, SINGLE PRECISION ARC-TANGENT FUNCTION
C00012 00006		[CONTINUATION OF ATAN]
C00014 00007	ATAN2 - FLOATING POINT, SINGLE PRECISION ARC-TANGENT WITH TWO ARGUMENTS
C00016 00008		[CONTINUATION OF ATAN2]
C00017 00009	SNCOS - F.P., SINGLE PRECISION SINE/COSINE USING TABLE LOOKUP
C00020 00010		[CONTINUATION OF EXECUTABLE CODE]
C00023 00011	   SINE AND COSINE LOOK-UP TABLES
C00027 00012		[CONT. SINE AND COSINE TABLES]
C00029 ENDMK
C⊗;
.TITLE ARITH

;Floating point square root.  Timing:  126 -- 135 microseconds.
;AC0 ← SQRT(AC0).  Clobbers AC1, AC2
;AC0:	argument
;AC1:	guess of answer
;AC2:	temp arith register

SQRTF:	ABSF 	AC0		;AC0 ← absolute value
	MOV   	R0,-(SP)
	STEXP 	AC0,R0		;R0 ← Exponent
	INC 	R0		 
	ASR 	R0		;R0 ← (Exponent +1)/ 2
	BCS 	SQR1		;If fraction ≥ half then branch.
	LDF 	#00234,AC1	;FRACTION[guess] ← .470 = sqrt(3/2)
	BR 	SQR2		 
SQR1:	LDF 	#00335,AC1	;FRACTION[guess] ← .673 = sqrt(3)
SQR2:	LDEXP 	R0,AC1		;EXP[Guess] ← exponent / 2
	MOV 	#3,R0		;R0 ← Number of Newton iterations
SQR3:	STF 	AC0,AC2		;temp ← arg
	DIVF 	AC1,AC2		;temp ← arg /guess
	ADDF 	AC1,AC2		;temp ← (arg/guess)+guess
	MULF 	#40000,AC2	;temp ← ((arg/guess)+guess) * .5
	STF 	AC2,AC1		;new guess ← temp
	SOB 	R0,SQR3		;Repeat iteration a few times
	STF 	AC1,AC0		;AC0 ← answer
	MOV	(SP)+,R0
	RTS 	PC		;Done
;ASIN - FLOATING POINT, SINGLE PRECISION ARC-SINE FUNCTION

;COMPUTES THE ARCSINE OF AN VALUE USING THE FOLLOWING ALGORITHM
;
;	ASIN(X) = ATAN(X/SQRT(1-X**2))
;
;WHERE X IS IN THE RANGE -1.0 TO +1.0.  THE ARGUMENT X MUST BE 
;LOADED INTO AC0 BEFORE THE CALL TO ASIN.  AFTER EXECUTION ASIN
;RETURNS THE ARCSINE VALUE IN AC0.  ASIN IS CALLED USING THE PC.
;RUN TIME IS APPROXIMATELY 380 MICRO SECONDS.

;REGISTERS USED:
;
;	AC0 PASSES ARGUMENTS
;	AC1,AC2,AC3  GARBAGED

;DEFINITIONS

ONE:	.WORD	40200,0		;FLOATING POINT 1.0
PIOT:	.WORD	40311,7733	;Pi/2

;START OF EXECUTABLE CODE

ASIN:	STF	AC0,AC3		;SAVE VALUE
	MULF	AC3,AC0		;GET X**2
	NEGF	AC0		;NOW -X**2
	ADDF	ONE,AC0		;NOW HAVE 1-X**2
	CFCC			;GET CONDITIONAL CODES
	BLE	ASN1 		;JUMP IF X WAS + OR -1.0
	JSR	PC,SQRTF	;COMPUTE SQUARE ROOT OF 1-X**2
	DIVF	AC0,AC3		;AC3= X/(1-X**2)**-.5
	STF	AC3,AC0
	JSR	PC,ATAN		;COMPUTE THE ARCTAN OF AC3
	RTS	PC		;RETURN
ASN1:	LDF	PIOT,AC0	;ANSWER WILL BE Pi/2 OR -Pi/2 DEG
	TSTF	AC3		;TEST SIGN OF ORIGINAL ARGUMENT
	CFCC
	BGE	.+4		;SKIP IF ARGUMENT POSITIVE
	NEGF	AC0		;ELSE SET ANSWER TO -Pi/2 DEG
	RTS	PC		;RETURN

;END OF "ASIN"
;ACOS - FLOATING POINT, SINGLE PRECISION ARC-COSINE FUNCTION

;COMPUTES THE ARC-COSINE OF AN VALUE USING THE FOLLOWING ALGORITHM
;
;	ACOS(X) = (PI/2) - ASIN(X)
;
;WHERE X IS IN THE RANGE -1.0 TO +1.0.  THE ARGUMENT X MUST BE 
;LOADED INTO AC0 BEFORE THE CALL TO ACOS.  AFTER EXECUTION ACOS
;RETURNS THE ARC-COSINE VALUE IN AC0.  ACOS IS CALLED USING THE PC.
;RUN TIME IS APPROXIMATELY 390 MICRO SECONDS.

;REGISTERS USED:
;
;	AC0 PASSES ARGUMENTS
;	AC1,AC2,AC3  GARBAGED


;START OF EXECUTABLE CODE

ACOS:	JSR	PC,ASIN		;COMPUTE ARCSINE OF ARGUMENT
	NEGF	AC0		;GET NEGATIVE OF ANSWER
	ADDF	PIOT,AC0	;NOW HAVE Pi/2 - ASIN(X)
	RTS	PC		;RETURN

;END OF "ACOS"
;ATAN - FLOATING POINT, SINGLE PRECISION ARC-TANGENT FUNCTION

;COMPUTES THE ARC-TANGENT OF A VALUE USING THE FOLLOWING ALGORITHM
;
;	ATAN(X) = X(B0+A1(Z+B1-A2(Z+B2-A3(Z+B3)**-1)**-1)**-1)
;	  WHERE  Z=X**2    FOR 0<X<=1
;		 Z=1/X**2  FOR 1<X
;
;THE ARGUMENT X MUST BE LOADED INTO AC0 BEFORE THE CALL TO ATAN.
;AFTER EXECUTION, ATAN RETURNS THE ARC-TANGENT VALUE IN AC0.  ATAN
;IS CALLED USING THE PC.  RUN TIME IS APPROXIMATELY 210 MICRO 
;SECONDS.  BEWARE!!!!!  BOTH + AND -∞ RETURN A VALUE OF PI/2.

;REGISTERS USED:
;
;	AC0 PASSES ARGUMENTS
;	AC1,AC2,AC3 GARBAGED


;DEFINITIONS:

A1:	.WORD	5046,23514	;2**-33
A2:	.WORD	73705,33431	;2**33
KB0:	.WORD	37462,154340	;0.1746554388
KB1:	.WORD	40730,61562	;6.762139240
KB2:	.WORD	40524,37327	;3.316335425
KB3:	.WORD	40271,66302	;1.448631538
KA1:	.WORD	40555,62164	;3.709256262
KA2:	.WORD	140743,65224	;-7.106760045
KA3:	.WORD	137607,107701	;-0.2647686202
TSIGN:	.WORD	0,0		;TEMP STORAGE FOR SIGN OF ARGUMENT

;START OF EXECUTABLE CODE

ATAN:	STF	AC0,AC1		;COPY ARGUMENT
	ABSF	AC1		;GET ABSOLUTE VALUE
	CMPF	A1,AC1		;COMPARE TO 2**-33
	CFCC			;TRANSFER CONDITIONAL CODES
	BLE	.+4		;SKIP IS X GREATER THAN 2**-33
	RTS	PC		;ELSE RETURN X AS ANSWER
	STF	AC0,TSIGN	;SAVE THE SIGN OF THE ARGUMENT
	CMPF	A2,AC1		;COMPARE X TO 2**33
	CFCC
	BGT	INRAN		;BRANCH IF X LESS THAN 2**33
	LDF	PIOT,AC0	;ELSE ATAN(X)= PI/2
	RTS	PC		;RETURN
;	[CONTINUATION OF ATAN]

INRAN:	LDF	ONE,AC0		;LOAD A F.P 1.0
	CLRB	TSIGN  		;USE TSIGN   TO INDICATE IF ABS(X).GT .1.0
	CMPF	AC1,AC0		;COMPARE ABS(X) TO 1.0
	MOV	R0,-(SP)	;SAVE REGISTER
	CFCC
	BLE	LE1		;BRANCH IF ABS(X) LESS THAN OR EQUAL TO 1.0
	DIVF	AC1,AC0		;ELSE REPLACE ABS(X) BY 1.0/ABS(X)
	MOV	#-1,R0		;LOAD ALL 1'S INTO R0
	XOR	R0,TSIGN	;COMPLEMENT SIGN AND SET .GT. 1 INDICATOR
	STF	AC0,AC1		;MOVE 1.0/ABS(X) BACK INTO AC1
LE1:	STF	AC1,AC3		;SAVE ARGUMENT
	MULF	AC3,AC1		;COMPUTE X**2=Z
	MOV	(SP)+,R0	;RESTORE REGISTER
	LDF	KB3,AC2		;PICK UP FIRST CONSTANT
	ADDF	AC1,AC2		;+ Z
	LDF	KA3,AC0		;PICK UP KA3
	DIVF	AC2,AC0		;NOW HAVE -A3/(B3+Z)   
	ADDF	AC1,AC0		;+Z
	ADDF	KB2,AC0		;NOW HAVE Z+B2 -(A3/(B3+Z))
	LDF	KA2,AC2		;PICK UP KA2
	DIVF	AC0,AC2		;DIVID -A2 BY PARTIAL SUM
	ADDF	AC1,AC2		;+Z
	ADDF	KB1,AC2		;+ B1
	LDF	KA1,AC0		;PICK UP A1
	DIVF	AC2,AC0		;DIVID A1 BY PARTIAL SUM
	ADDF	KB0,AC0		;+B0
	MULF	AC3,AC0		;MULTIPLY BY ORGINAL ARGUMENT
	TSTB	TSIGN   	;CHECK IF ABS(X) WAS GREATER THAN 1.0
	BEQ	.+6		;SKIP IF ABS(X) LESS THAN 1.0
	SUBF	PIOT,AC0	;ELSE ATAN(X)=-(ATAN(1/X)-PI/2)
	TST	TSIGN		;CHECK SIGN OF X
	BPL	.+4		;SKIP IF (SIGN.XOR.GT1) IS POSITIVE
	NEGF	AC0		;ELSE ATAN(X)=-ATAN(X)
	RTS	PC		;RETURN

;END OF "ATAN"
;ATAN2 - FLOATING POINT, SINGLE PRECISION ARC-TANGENT WITH TWO ARGUMENTS

;COMPUTES THE ARC-TANGENT OF A/B USING THE ATAN ROUTINE.  HOWEVER, SINCE
;TWO ARGUMENTS ARE USED SINGULARITIES ARE AVOIDED AT MULTIPLES OF PI/2 
;AND THERE IS NO AMBIGUITY CONCERNING QUADRANTS.  THE ARGUMENT A MUST
;BE LOADED INTO AC0 AND B INTO AC1 BEFORE CALLING ATAN2.  AFTER EXECUTION,
;ATAN RETURNS THE ARC-TANGENT VALUE IN AC0.  RUN TIME IS APPROXIMATELY
;250 MICRO SECONDS.

;REGISTERS USED:
;
;	AC0,AC1 PASSES ARGUMENTS
;	AC2,AC3 GARBAGED

;DEFINITIONS

TSIGN2:	.WORD	0,0
PI:	.WORD	40511,7733	;Pi

;START OF EXECUTABLE CODE

ATAN2:	STF	AC0,AC2		;GET ABSOLUTE VALUE OF A
	ABSF	AC2
	STF	AC1,AC3		;GET ABSOLUTE VALUE OF B
	ABSF	AC3
	CMPF	AC2,AC3		;COMPARE ABS(A) TO ABS(B)
	CFCC
	BPL	REVAB		;BRANCH IF ABS(A)> ABS(B)
	STF	AC1,TSIGN2	;ELSE SAVE SIGN OF B FOR LATER
	DIVF	AC1,AC0		;COMPUTE A/B
	JSR	PC,ATAN		;COMPUTE TAN(A/B)
	TST	TSIGN2		;CHECK THE SIGN OF B
	BPL	ATDNE		;IF B>0, LEAVE SGN(ATAN)=SGN(A) AND EXIT
	TSTF	AC0		;TEST IF A POSITIVE
	CFCC
	BGT	ATAN2A		;BRANCH IS A POSITIVE
	ADDF	PI,AC0		;NO, SECOND QUADRANT, ADD PI
ATDNE:	RTS	PC		;RETURN
ATAN2A:	SUBF	PI,AC0		;YES,3RD QUADRANT, SUBTRACT PI
	RTS	PC
;	[CONTINUATION OF ATAN2]

REVAB:	NEGF	AC1		;IF ABS(A)<ABS(B) THEN COMPUTE -B/A
	STF	AC0,TSIGN2	;SAVE THE SIGN OF A
      	DIVF	AC0,AC1
	STF	AC1,AC0		;LEAVE -B/A FOR ATAN
	JSR	PC,ATAN		;COMPUTE ATAN(-B/A)
	TST	TSIGN2		;CHECK THE SIGN OF A
	BMI	XYZ		;BRANCH IF A NEGATIVE
	ADDF	PIOT,AC0		;ELSE ADD TWO PI TO ATAN
	RTS	PC		;RETURN
XYZ:	SUBF	PIOT,AC0		
	RTS	PC

;END OF "ATAN2"
;SNCOS - F.P., SINGLE PRECISION SINE/COSINE USING TABLE LOOKUP

;THIS PROGRAM  CALCULATES  BOTH THE  SINE AND  THE COSINE  OF A  GIVEN
;ANGLE. THE  RESULTING VALUES ARE ACCURATE  TO WITH IN +  OR - .0003 .
;RUN  TIME IS  APPROXIMATELY  142  MICROSECONDS  PER  CALL.    ADD  40
;MICROSECONDS IF THE ANGLE IS GREATER THAN + OR - 4Pi. THIS ROUTINE IS
;SOMEWHAT  FASTER  THAN  THE  SYSTEM  SINE,  COSINE  EXPANSION FORMULA
;ROUTINES SINCE A TABLE LOOK-UP TECHNIQUE IS EMPLOYED.   THERE ARE TWO
;ENTRY POINTS  INTO THIS ROUTINE: SNCOSR  AND SNCOSD.   SNCOSR IS USED
;FOR ANGLES GIVEN IN RADIANS AND  SNCOSD FOR ANGLES IN DEGREES.   BOTH
;ROUTINES ARE CALLED USING THE PC.  THE ANGLE  TO BE CONVERTED MUST BE
;LOADED INTO AC0  BEFORE THE CALL TO SNCOS.  AFTER EXECUTION, THE SINE
;OF THE ANGLE IS LEFT IN AC0 AND THE COSINE IN AC1.

;REGISTERS USED:
;
;	AC0,AC1 PASS ARGUMENTS
;	AC2,AC3 GARBAGED

;DEFINITIONS

SNR==%2
CSR==%3
DEGCV:	.WORD	41466,5541	;'40000/360.0
RADCV:	.WORD   43042,174604	;'40000/2*Pi
SCNORM:	.WORD   34600,00000	;1/'40000
SCINT:	.WORD	43600,00000	;'40000
SIGN:	0

;START OF EXECUTABLE CODE

SNCOSD:	MULF	DEGCV,AC0	;NORMALIZE, 0 TO 360 = 0 TO 40000
	BR	SVREG

SNCOSR:	MULF	RADCV,AC0	;NORMALIZE, 0 TO 2Pi = 0 TO 40000
SVREG:	MOV	R0,-(SP)	;SAVE REGISTERS
      	MOV	R2,-(SP)	
	MOV	R3,-(SP)
ANNORM:	STCFI	AC0,SIGN	;GET THE NORMALIZED ANGLE
	CFCC			;TRANSFER THE CONDITION CODES
	BCC	RTRGE		;SKIP IF NO CONVERSION OVERFLOW
	MODF	SCNORM,AC0	;NORMALIZE BETWEEN +0 AND 1
	MULF	SCINT,AC0	;SCALE TO 0 TO 40000
        BR 	ANNORM		;GO TRY INTEGERIZING AGAIN
	[CONTINUATION OF EXECUTABLE CODE]

RTRGE:	MOV	SIGN,R0		;LOAD INTEGER
	BGE	.+4		;GET ABSOLUTE VALUE
	NEG	R0
	ASH	#2,R0		;MULT BY 4 TO USE AS F.P. INDEX
	MOV	R0,R3		;NEED THIS LATER
	BIC	#177400,R0	;6 LSB BECOME FINE ESTIMATE OF SIN, COS
	LDF	COSF(R0),AC0	;FINE ESTIMATE OF COS TO AC0
	SWAB	R3		;POSITION NEXT 6 BITS TO USE AS INDEX
	ASH	#2,R3		;x 4 FOR ANOTHER F.P. INDEX
	MOV	R3,R2		;SAVE REMAINING BITS
        STF	AC0,AC1		;FINE ESTIMATE OF COS TO AC1
	BIC	#177400,R3	;NOW HAVE INDEX INTO ROUGH ESTIMATE OF S,C
    	BIT	#400,R2		;TEST IF THETA IN QUADRANT 2 OR 4
	BNE	Q2OR4		;BRANCH IF IT IS
	LDF	RUF(R3),SNR	;ROUGH ESTIMATE OF SINE TO SNR
	NEG	R3		;GET INDEX FROM OTHER END OF LIST
	LDF	RUFE(R3),CSR	;  "      "      " COSINE TO CSR
        BR 	CHKS
Q2OR4:	LDF	RUF(R3),CSR	;REVERSE VALUES IF IN QUAD 2 OR 4
	NEG	R3
	LDF	RUFE(R3),SNR
CHKS:	BIT	#1000,R2	;COMPLEMENT SIN-R IF IN QUAD 3 OR 4
	BEQ	.+4
	NEGF	SNR
	MULF	SNR,AC0		;COMPUTE SIN/COS PRODUCTS  AC0=SNR*CSF
	SWAB	R2		;CHECK IF IN QUADRANT 2 OR 3
	INC	R2
	MULF	SINF(R0),SNR	;  SNR=SNR*SINF
	BIT	#2,R2		
	BEQ	NOC		;SKIP IF THETA IN QUAD 1 OR 4
	NEGF	CSR		;ELSE COMPLEMENT ROUGH COSINE
NOC:	MULF	CSR,AC1		;  AC1=CSR*CSF
	MOV	(SP)+,R3	;RESTORE REGISTER
	MULF	SINF(R0),CSR	;  CSR=CSR*SNF
	MOV	(SP)+,R2	;RESTORE REGISTER
	ADDF	CSR,AC0		;SIN THETA = SNR*CSF + CSR*SNF
	MOV	(SP)+,R0	;RESTORE REGISTER
	TST	SIGN		;COMPLEMENT SIN IF THETA NEGATIVE
	BGE	.+4
	NEGF	AC0
	SUBF	SNR,AC1		;COS THETA = CSR*CSF - SNR*SNF
	RTS	PC		;EXIT

;END OF "SNCOS" EXECUTABLE CODE
;   SINE AND COSINE LOOK-UP TABLES

RUF:

	.WORD       0,     0, 36711,  5260, 37110,175460, 37226,124405
	.WORD   37310,136466, 37372,131163, 37426, 40204, 37457, 10242
	.WORD   37507,142702, 37540, 56024, 37570,147714, 37610,107223
	.WORD   37624,120062, 37640,115345, 37654, 76324, 37670, 42052
	.WORD   37703,167425, 37717, 75712, 37732,164200, 37746, 31565
	.WORD   37761, 55352, 37774, 56447, 40003,116075, 40010,172633
	.WORD   40016, 34732, 40023, 64052, 40030, 77700, 40035, 77721
	.WORD   40042, 63631, 40047, 33126, 40053,165512, 40060,102673
	.WORD   40065,  2363, 40071, 64102, 40075,127371, 40101,154161
	.WORD   40105,162004, 40111,150423, 40115,117403, 40121, 46475
	.WORD   40124,155462, 40130, 44123, 40133,112032, 40136,137006
	.WORD   40141,142630, 40144,125131, 40147, 65730, 40152,  4647
	.WORD   40154,101537, 40156,154236, 40161,  4410, 40163, 12107
	.WORD   40164,175013, 40166,135007, 40170, 51770, 40171,143635
	.WORD   40173, 12276, 40174, 35450, 40175, 35254, 40176, 11443
	.WORD   40176,142155, 40177, 47155, 40177,130417, 40177,166103
RUFE:	.WORD   40200,     0

SINF:

	.WORD       0,     0, 35311,  7733, 35511,  7733, 35626,145744
	.WORD   35711,  7733, 35773, 51722, 36026,145733, 36057,166722
	.WORD   36111,  7706, 36142, 30671, 36173, 51651, 36212, 35314
	.WORD   36226,145701, 36243, 56265, 36257,166650, 36274, 77231
	.WORD   36311,  7610, 36325,120165, 36342, 30540, 36356,141111
	.WORD   36373, 51460, 36403,171012, 36412, 35173, 36420,101352
	.WORD   36426,145530, 36435, 11705, 36443, 56060, 36451,122232
	.WORD   36457,166402, 36466, 32550, 36474, 76715, 36502,143057
	.WORD   36511,  7220, 36517, 53357, 36525,117514, 36533,163647
	.WORD   36542, 27777, 36550, 74126, 36556,140252, 36565,  4374
	.WORD   36573, 50514, 36600,146315, 36603,170362, 36607, 12426
	.WORD   36612, 34471, 36615, 56532, 36620,100572, 36623,122631
	.WORD   36626,144666, 36631,166722, 36635, 10754, 36640, 33005
	.WORD   36643, 55034, 36646, 77062, 36651,121106, 36654,143130
	.WORD   36657,165151, 36663,  7170, 36666, 31205, 36671, 53221
	.WORD   36674, 75232, 36677,117242, 36702,141250, 36705,163255
;	[CONT. SINE AND COSINE TABLES]

COSF:

	.WORD   40200,     0,40177,177777,40177,177773,40177,177765
	.WORD   40177,177754,40177,177741,40177,177724,40177,177704
	.WORD   40177,177661,40177,177634,40177,177605,40177,177553
	.WORD   40177,177516,40177,177457,40177,177416,40177,177353
	.WORD   40177,177304,40177,177233,40177,177160,40177,177103
	.WORD   40177,177023,40177,176740,40177,176653,40177,176563
	.WORD   40177,176471,40177,176375,40177,176276,40177,176175
	.WORD   40177,176071,40177,175762,40177,175652,40177,175537
	.WORD   40177,175421,40177,175301,40177,175156,40177,175031
	.WORD   40177,174701,40177,174547,40177,174412,40177,174253
	.WORD   40177,174112,40177,173746,40177,173600,40177,173427
	.WORD   40177,173254,40177,173076,40177,172715,40177,172533
	.WORD   40177,172345,40177,172156,40177,171764,40177,171567
	.WORD   40177,171370,40177,171167,40177,170763,40177,170554
	.WORD   40177,170343,40177,170130,40177,167712,40177,167472
	.WORD   40177,167247,40177,167021,40177,166572,40177,166340